#include "uhm.hxx"

#define N 1

using namespace uhm;

typedef ID_Node_H_<N>          Id_Nd;
typedef ID_Element_<N>         Id_El;

#ifdef NOPIV
typedef CHOL_nopiv_Linal_        Uhm;
#endif

#ifdef PIV
typedef CHOL_piv_Linal_          Uhm;
#endif

#ifdef INCPIV
typedef CHOL_incpiv_Linal_       Uhm;
#endif

typedef Direct_Solver_<Id_Nd,Id_El,Uhm> Solver;

int main (int argc, char **argv) {

  initialize(argc, argv);

  // ** Environments
  Int_
    type = Mat_Base_::DOUBLE,
    uplo = Mat_Base_::LOWER_TRIANGULAR,
    sym  = Mat_Base_::SYM,
    rhs  = 1,
    blk  = 256;

  Int_
    kind = Mesh_Base_::NODE_DEFAULT;


  // ** Solver initialization
  Solver solver(type, uplo, sym, rhs, blk);

  // ** Import a mesh from a file_in
  option_begin();
  std::string fin;
  get_option_string("-in", "Reading a file", __FILE__,
                    "../../../../uhmfile/toy_1.uhm", fin);
  Double_ err;
  get_option_double("-err", "Error level to determine PASS/FAIL", __FILE__,
                    1.0e-3, err);
  Int_ caching;
  get_option_int("-caching", "Setting caching level", __FILE__,
                 6, caching);
  solver.set_caching(caching);

  Int_ n_iter;
  get_option_int("-iter", "Number of solving iteration", __FILE__,
                 2, n_iter);

  option_end();

  Int_ r_val;
	
	  time_begin();

  std::cout << argv[0] << "  Opening file ... " << fin << std::endl;
  std::cout << "  Caching ... " << caching << std::endl;
  
  std::ifstream in;

  in.open(fin.c_str());
  stream2mesh<Solver::_Mesh_>(in, solver.mesh());
  in.close();
  
  // ** Create random matrices in work
  UHM_ERROR((mesh2cache<Solver::_Mesh_,Solver::_Cache_>(solver.mesh(), solver.work())),
            ">> Error in mesh2cache");
  
  solver.work().set_matrices(type, rhs, blk);
  solver.work().create_matrices();
  solver.work().random_matrices(true);

  // ** Clear connectivity for successive iteration
  UHM_ERROR(solver.clear(), ">> Error in Clearing");
  
  for (Int_ i=0;i<n_iter;++i) {

    in.open(fin.c_str());
    stream2mesh<Solver::_Mesh_>(in, solver.mesh());
    in.close();

    // ** create tree and allocation
    time_in("direct:: Creating Solver (analysis)");
    UHM_ERROR(solver.create(),    ">> Error in Creating");
    time_out();

    time_in("direct:: Decomposing");
    UHM_ERROR(solver.updecompose(), ">> Error in Decomposing");
    double t = time_out();

    Double_ flop, recorded;
    solver.get_flop(flop, recorded);

    std::cout << "\n\nFLOP = " << (flop/1.0e9)
              << ", Updated  = " << (recorded/1.0e9) << "  GFLOP\n";
  
    std::cout << "\n\nFLOPS= " << (flop/1.0e9)/t
              << ", Modified = " << ((flop-recorded)/1.0e9)/t << "  GFLOPS\n";

    Double_ storage, max_storage;
    storage_counter(storage, max_storage);
    std::cout << "\n\nStorage = " << (storage/1.0e6) 
              << ", " << (max_storage/1.0e6) << "  MBytes\n\n";


    time_in("direct:: Solving");
    UHM_ERROR(solver.solve(), ">> Error in Solving");
    time_out();

    time_in("direct:: Checking");
    Double_ res;
    UHM_ERROR(solver.check(res), ">> Error in Checking");
    time_out();

    std::cout << "\n\nResidual = " << res << "\n\n";
    r_val += (res < err ? 0: -1);
    std::cout << (res < err ? "[PASS]\n\n" : "[FAIL]\n\n") ;

    UHM_ERROR(solver.clear(), ">> Error in Clearing ");
    std::cout << "  " << i << "  iteration ....\n";
  }
  
  UHM_ERROR(solver.flush(), ">> Error in Flushing");
  
  time_end();
  
  finalize();

  r_val = (r_val < 0 ? -1 : 0);

  return r_val;
}


